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SYSTEM EVALUATION FROM IMPULSE RESPONSE FRAGMENTS 


INTRODUCTION 


Estimation or identification of transfer functions is often required 
in control techniques or any other field in which the knowledge of system dy- 
namics is a necessity. In aerospace techniques, it might be desirable to 
identify the flight dynamics of a vehicle. The two main cases of application 
are (1) identification of vehicle dynamics in self-adapting or self-learning 
control procedures and (2) confirmation of vehicle dynamics in flight testing 
of newly developed aerospace vehicles. 

Different approaches of estimation and identification have been made, 
but all of them are based on statistical methods. A whole family of procedures 
ends up with the impulse response of a system to be identified. Because of 
the inherent characteristic of these procedures, the im^plse response is 
represented by an equally sampled time series. The description of the sys- 
tem by its impulse response is often not desirable for data handling. In many 
cases, e.g. , in control law optimization, it is more convenient to have the 
vehicle dynamic described by a rational expression of polynominals [ 1] . 

Furthermore, the impulse response represented by the time series 
is often available only for a limited time interval. This can occur if the 
system under investigation is either instable, indifferent, or only slightly 
damped. In these cases the calculation of the impulse response has to be 
stopped after a time interval to be defined. Another reason for only a limited 
knowledge of the impulse response might be a limited computer capacity 
available for its calculation. 

Given a sampled record of a portion of the impulse response of a lin- 
ear system, estimate the transfer functions, G (S) and G (z), of the system. 

A procedure which covers the preceding requirements is found by use of 
rational approximation. The impulse response described by a time series is 
converted to a rational expression in z. After determination of the poles and 
zeros of the function, a transformation to the s-domain can then be provided 
if required. In the following section, the appropriate procedure is described, 
and restrictions are discussed. A digital program for implementation of the 
method is attached. 
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TASK DESCRIPTION 


By use of some identification or estimation procedure based on 
statistical theory, up-to-then totally unknown systems have to be identified 
by its impulse response. The impulse response is represented by a time 
series equally sampled by a sampling interval T, which is known from the 
identification procedure (Fig. 1). As a shortcoming, the impulse response is 
known only over a limited time interval. 



Figure 1. Time series represen- 
tation of impulse response. 


The available data must be 
converted to a form which describes 
the vehicle dynamic completely in 
the sampling moments and which 
can be used for further data pro- 
cessing. 

THE RATIONAL 
APPROXIMATION 

The z -transform of the time 
series x(kT) is defined to be 

a - k 

G(z) = c +c z~ + . . . + c, z +... 


( 1 ) 


where c = x(kT) with sampling interval T. 

K 


The desired equivalent rational expression in z is of the form 


F(z) 


a o + a i z_1 + 


. . + a z 
m 


-m 


V 


b^z 1 + 


b z 
n 


-n 


(2) 


If (1) is any series in z , there exists a unique rational "approxima- 
tion" of the form (2) for each pair of integers m and n that agrees when ex- 
panded term by term with the series for more terms than any other rational 
expression with smaller or equal m and n [2] . It follows that, if G(z) has a 
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rational expression, as is the case with most engineering systems, the exact 
rational expression should be found when m and n are chosen large enough. 

The choice of m and n does create a problem which will be discussed later. 

Once m and n are chosen, the rational approximation is found in multi- 
plying G(z) by the denominator of F(z) and collecting the like powers of z“ 


-1 -2 -k 

c. + c,z + cz + . . . + c, z +... 
0 12 k 


, ,-1,-2 , -n 

b n + b,z + b.z + . . . + b z 
0 12 n 


Vo + Vi 2 " 1 + b oV ‘ 2 + • • • + b oV" k 


, -1 , -2 . -3 , -(k+1) 

^1 C Q Z + ^1 C 1 Z + ^1 C 2 Z + • • • + 


(3) 


Vo + ( Vl + b l C 0 )z_1 + ( V 2 + Vl + Vo )Z ’ 2 

The first m+1 powers of z 1 can be forced to agree with the numerator 
of F(z) , and the next n can be forced to vanish. 

The equations for this are 


a k 


E cb - 

H , ii 


k ^ m 


(4) 


i+j=k 


thus, 


a o = Vo 


a ! ’ Vi + Vo ' 


(5) 
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and 


y. b.c. = 0 

1 — 1 l 1 

i+j=k , k = m+1, m+2, ... m+n 

i<n 

with n 2: m; 
thus, 


b A c ^ + b.c + . . . + b ,c. = 0 
0 m+1 1 m m+1 0 


be + b c + ...+b _c n = 0 
0 m+2 1 m+1 m+2 0 


( 6 ) 


b A c + b,c , + ...+ b c. = 0 
On 1 n-1 n 0 


(7) 


be _ + be + . . . + b c = 0 
0 n+1 In n 1 


b.c + b.c , + ... be 
0 n+m 1 n+m-1 n m 


0 


The last equations are solved for b., b,, ... b and then substituted in the 

0 1 n 

first m+1 to obtain a A , a,, ... a . Clearly, one of the b. can be chosen 

0’ I’ m J l 

arbitrarily (anything except 0) , and the rest are unique if the system of equa- 
tions is not singular. If m and n are chosen too large, the system of 
equations will be singular. 

If the rank of the system is r, the denominator of F(z) should be of 

degree r in z *. If the system is solved using n too large, theoretically the 
numerator and denominator of F(z) will have common factors which could 
be factored out. In practice, these factors may not be exactly the same be- 
cause of the roundoff errors. 
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Some procedure for checking whether the value of n is too large should 
be used to prevent the introduction of zeros in the denominator and also to 
prevent the solution of a singular system of equations. 


SCALING 

Assuming that the time sequence x(nAt) has been received from the 
time function x(t) through an ideal sampler, x(nAt) is attenuated by a 

factor — [3]. For correct scaling, the numerator of f(z) must be multiplied 

by At. 


RECONVERSION TO THE TIME DOMAIN 


For later use in the realization of the procedure described previously, 
it is presumed [4] that F(z) is reconverted to the time domain by 

y. = a„x, + a„x, _+...+ a x, _ - b,y, _ . . . -b y, ^ . . 

Ok 1 k-T m k-mT rk-T nTk-mT (8) 

where x and y are the input and output time sequences of the digital filter, 

K K 

and x^ - mT and y^ - mT are the m- T earlier values of the input and output 

time sequences. Applying an impulse as an input to equation (8), the impulse 
response of F(z) can be calculated as a time series representation. 

IMPLEMENTATION OF THE PROCEDURE IN A DIGITAL COMPUTER 


For implementation of the method just described, it has to be realized 
that, according to the task description, the vehicle dynamics are completely 
unknown; thus, information on the order of the rational expression in both 
numerator and denominator's not available. Consequently, a scheme must 
be established to evaluate m and n. A suitable scheme has been found by 
taking logical steps as follows: 
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1. The order of the numerator and denominator of equation (2) is 
assumed to be equal; thus m = n. As later shown, this assumption does not 
introduce a significant error. 

2. An initial estimate on m is made assuring that 
m estimate ~ m actual 

3. With this estimate in m, the procedure previously described 
will be executed. If the estimate was too high, the matrix which has to-be 
solved in the process of the procedure will be ill conditioned or will provide 
an otherwise insignificant solution. 

4. To confirm the result received from step 3, the coefficients of 
F(z) will be checked out by calculating its impulse response according to 
equation (8). The n+m+1 values of the impulse response time series will be 
calculated. If these values disagree with the input time series, m is reduced- 
by 1 and the complete procedure is repeated until correspondence between 
both time series is received. 

5. If we assume that the sampling rate has been chosen at least 
four times higher than the highest natural mode of the system, all roots of 
F(z) will be located in the right half of the z plane. This will result in alter- 
native signs for the polynomial of F(z). Thus, an alternative sign condition 
might be established and introduced prior to the implementation of step 4. 
However, this condition is not necessarily acceptable or useful in all cases 
of application. 

This scheme is, of course, sensitive to the accuracy of the calcula- 
tion. Thus, the number of digits which have to agree has to be evaluated 
taking into account the number of significant digits through the whole sequence 
of calculations. Furthermore, in following this procedure, possible dis- 
turbances by noise of the initial input data have to be considered. 

Additional confidence in this described checkout procedure may be 
established by extending the number of values from the time series to be 
compared, e. g. , 2 • (n+m+1). 

A digital program performing the procedure described is presented 
in the appendix. 


6 



APPLICATION OF THE DIGITAL PROGRAM 


As an example and to show the significance of the method, the digital 
program is applied to an impulse response of which the transfer function 
in s is known. To provide exact knowledge to the input data, the impulse 
response has been received by inverse Laplace transformation of the transfer 
function in s. 

The example has been chosen to be 


(1 + 


F 1 (s) = 


31.4 


s) 


, . 2 12 , 

(l + r-7T s + s ) 


6.28 


6. 28 


(9) 


with co = 6. 28 and £ = 1. 

The impulse response is printed in Table 1 for t = 0. 01 sec. Figure 2 
shows the impulse response as a curve. 

From the given example, we know that m = 2. Assuming that this is 
not known, we make an initial estimate of m = 4; thus, the first m + n + 1=9 
values are used as an input to the program. Futhermore, m = n and t = 0. 01 
sec. 


Executing the program results in 

„ , , 0. 1256 • 10 -1 - 0. 882432663955 • 10 -2 z _1 + 0. 133715714656 • 10~ 5 z~ 2 

F (z) = 

1 0. 1 • 10 -1 - 0. 187758651589 • 10 1 z _1 0.881226462775 • 10° z -2 

with gain factor GF = 1. 02666622147 . 

Table 2 shows the impulse response derived from the calculated F^(z) 
according to equation (8.). Comparing with the input response of Table 1 
and Figure 2, it can be seen that a slight deviation builds up with increasing 
numbers of time sequence values, due to roundoff errors, but that it still 
does not exceed 15 percent of the hundredth already very small input value. 
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TABLE 1. IMPULSE RESPONSE VALUES AND ROOTS FOR EXAMPLE F , (s) 
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Figure 2. Co 


TABLE 2. IMPULSE RESPONSE TIME SEQUENCE, RECONVERTED FROM FjU). 
POLES AND ZEROS OF F 1 (z) AND CORRESPONDING ROOTS IN s-DOMAIN 
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RELATION BETWEEN z- AND s-TRANSEER FUNCTION 


Considering the theoretical approach of the procedure and its execution 

as previously described, it can be stated that an exact z-transform has been 

provided with respect to F , . . The accuracy is limited only by the number 

Is / 

of significant digits used in the course of the calculation. 

Assuming that the system to be identified has been a digital one, it is 
totally described by the time series which has been used as input. Yet it has 
to be realized that, in most of all cases, specifically in identifying the dynamic 
of a vehicle, the system is an analog one. Thus, the dynamic of the calculated 
F^ differs from the actual F^ represented by the impulse response. 

For further consideration, we assume that the sampling rate has been 
chosen as high and that with respect to the natural frequency of the system, 
its dynamical behavior between sampling instances can be neglected. 

Any remaining difference in the dynamical behavior between the calcu- 
lated F^ and the actual F ^ is largely dependant on the sampling rate. 

If the sampling rate has been selected sufficiently high with respect to the 
natural modes of the system and the interesting frequency range, these dif- 
ferences might be negligible. 

To cover more demanding requirements, a transformation to the 
s -domain has to be provided. 

Thus, after having solved F(z) for poles and zeros, the appropriate 
roots in s are found according to Figure 3 as 


= *\|(Re z)^ + (lmz)^ 


Re s 


A. 

T 


In z 
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Figure 3. Correspondence between z- and s-domain. 


and 


Im S 


tan 


-1 Im z 
Re z 


T 


The gain factor of the system is received out of F(z) by letting 

z *_ 1 . 


The executive of the transformation to the s-domain has been provided 
for the example previously discussed. The printout in Table 2 shows that 
some deviation from the original transfer function in s exists. This deviation, 
however, would not be significant in a practical case of application. Further- 
more, it can be seen that the assumption m = n does not result in any signifi- 
cant error. The erroneously introduced second zero at -1317 . . . has no 
practical effect on the dynamic of the system. 
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ACCURACY 


In working with the procedure outlined in this report, it must be con- 
sidered that it is sensitive 

1. To noise superimposed to the impulse response as a result of 
errors in the identification procedure, and 

2. To the accuracy of the execution in the digital program. 

Both effects become more effective with decreasing T because of the 
decreasing gradient between sample values. Thus, noise becomes more and 
more significant and may create necessity for the application of smoothing 
procedures. Furthermore, decreasing T requires increasing accuracy, i. e. , 
number of significant digits of the execution procedure. 

Inaccuracy effects become even more significant for more complex 
transfer functions where the number of operations, necessary for the execu- 
tion of the program, increases. 

CONCLUSION 


The procedure outlined in this report allows the exact z transformation 
from a continuous time function presented by an equally spaced time series 
and conversion of the resulting power series in z to a rational expression in 
z. By a proper choice of the sampling interval, the dynamical behavior of 
the transfer function in s can be significantly approximated by the transfer 
function in z with respect to practical application. If necessary, however, a 
transformation to the s-domain is easily provided. By introduction of a special 
program subroutine and by taking an initial estimate on the order of the trans- 
fer function in s, it is possible to apply the procedure even if the order of the 
system under investigation is not known. It has to be made sure only that the 
estimate is greater or at least equal to the actual value. For the implementa- 
tion of the procedure, only a limited number of values of the impulse response 
time series have to be known. In the ideal case, this number equals 1 plus 
twice the order of the denominator of the transfer function in s. For practical 
application, it might be useful to have a slightly higher number for checkout 
purposes. In any case, however, only a fraction of the complete impulse 
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response of the system under investigation has to be known. Thus, the pro- 
cedure is suitable for data processing in connection with the identification of 
unknown plants. 

It must be realized that the digital program in the Appendix has been 
developed to prove the practical feasibility of the described method. It is in 
no way optimal with respect to the required computer capacity. If the practical 
application of the described method is envisaged, it- is suggested that the digital 
program be reviewed. It is especially recommended to increase the accuracy 
of the procedure. Furthermore, it might be useful to extend the checkout 
procedure with respect to the number of time series values to be compared. 

In addition, the integration of a smoothing procedure for smoothing of the 
input values as a fixed subroutine of the program should be considered. 


George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, Alabama 35812, June 1, 1971 
70-103-19-04 
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APPENDIX 


DIGITAL PROGRAM 


IMPLICIT DOUbLE PRECISION (A-H.O -2) 

DIMENSION C C&OO ) * A (i»UO ) tB (SOO ) 

DIMENSION E ( bOU ) t D ( 1 02 ) 

DIMENSION p ( 800 ) 

DIMENSION NTK(SOU).NTIISOO)»T6I50U)»G(SOO)»^(bOO),Y(SCO> 
DIMENSION TOO t sot) r ,00 15001 ,XX(500l,YYt500) 

C RATIONAL APPROAIMATJON in DOUBLE PRECESION 
NAME LIST /InPuT/ M»N,DT,C 

1 READIS, INPUT) 

00103 

2 N=N- t 
M®h“ 1 

004 1*1,30 
B ( 1 ) *0 • UDO 

4 CONTINUE 
N 2 ■ N 2 * 2 

3 M N ■ M N 
M 1 *M* 1 
N I « N ♦ I 

N 2 ■ M ♦ M ♦ I 
rtKITL(fc,IS)M,N 

IS FORMAT ( / / 3 H M«jS,/,3h N»IS/) 
i r E l 6 , 123) OT 
123 pORMA|(4H UT«10F4/) 

IF IM«-I..Eo.Q> 00 10 bU 

CALL PRINT ( 12HINPUT DATA ,N2,C> 

c simultaneous equations are foKmeu 

M fv b h ♦ N 

K*D 

Mn»2 

5 LL® t MN/2 ) +MM 
002S J*l , N 1 

K = K ♦ 1 

E IK >»C ILL ) 

Ll»LL- l 
25 CUNT INUE 

1 P (MM ,NE ,H 1 ) 00 TO 9 
00 TO 8 
V M (i a M M ♦ 1 
0 0 TO 5 

8 CALL PRINT ( 12M E ,3U,E> 
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C C 0 £ F • OF 0 A fit: S I Ori£0 in £, aNO PUT IN THE FOKm REwUI«LO t) T oSPE 
C SUdHUUTlNE 

K-U 

l)U 10 J- 1 ,M 
L--M+ J 

0010 1 = l til 
K*l(+I 

L * L ♦ H * 1 

20 MiO=ElL> 

C = -M 

0030IJ J- 1 , M 
K -<*1 
L — L ♦ M + I 

3U0 F<K»=-E(L) 

CALL PRINT ( I 2h F ,30, F) 

k=m** 2 +m 

Call opse i y » n ,n ,det ) 
c a is changed into oou 3 L e precesion 

| 00 003UN I»1 |N 

j-M-N-N • ( I - 1 ) 

3 0 H a l I » «F ( J J 
D U H b I - l , N 
A * N - I ♦ l 

SS BlK+l)«SlK) 

a ( i »« i .u 

call printu2h a array ,ni,h) 
c a array is calculates 
a < i i = c < i » 

Ull >*1 .0 
ooaj 1-2, Ml 
11-1 
T-0.0 

0030 J- 1 , I 
r = T+d ( J I -C ( l I I 
30 11*11—1 

ia miui 

c impulse is calculaTeo as a Checkout for the program 

011 02) *0*0 
0 l 1 ) -A l 1 ) 

0 0999 1-2,101 

u- I 

T-0,0 

JO d dH J - 2 , N 1 

U-L- 1 

A-L 

iP<u.or,ii k*io2 
r--B i j) *uno *t 
aaa continue 

IFd.GT.Nn A ( 1 ) -U • o 
0 ( I ) -A 1 1 ) *T 
*99 CONTINUE 
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00777 I » 1 i N 1 
A I 1) >A ( I >*0T 
777 C u N T 1 Nli £ 

CALL PH INT < 12HA AHH Ar 

K a 0 • 0 

T»U.O 

001222 1 = 1 » t* 1 
T = T ♦ A ( I I 
R*R*B I l > 

1222 CON I l NUE 

of = r /n 

(iK 1 I E l 6 , 901 ) GF 
9 U 1 F 0 ri .1 A I ( 9 H 6 F“, 029.12/) 

TOL-0.0 

OUbbb 1=1,01.2 

I F < B 1 1 ) .LT . TQL .OH • b < 1 ♦ 1 > «GT • TUL ) 60 TO 2 
bbb CONTINUE 

call print ( 1 2H b stable .m.b) 

THE IMPULSE lb C HE C < E 0 MfH Th£ I M P UL SE -H E 5 P 0 N SE INPUT 
00111 1 = 1 .nZ 

I f ( ABS ( 0 < I 1 -C ( l I I • g£ • 0 . uoou I ) 60 TO 2 
III CONTINUE 

29 Call PRINT (I2H 0 IMPULSE ,101,01 

K * .V 

0 U 1 u I s 1 , N 
AUlsAIKI/AlII 
10 K*K-1 

CALL «TPOLT<N»A.bO,I.E-2b.RTR,RTl.CONV,A»B.C,O.E> 

a H 1 T E ( £ , 032 ) 

j 3 2 FORMAT I • NU« ROOTS Z») 

0 0 9 9 2 I » l , N 

AHl I E 1 6 . 9 9 I ) RTMII.KTIUI 
9 9 1 FUaMATIJH 2« 02H.I2.I1H ♦ AND - 0 02H.I2) 

h 9 2 continue 

UU9 93 1 = 1 ,N 

T O I 1 ) «DSORT IRTk m »RTR ( 1 )*rT 1 I I ) »RT H l I ) 

6 < I l=RT I ( I J /HTH ( I I 
X ( I ) =UAT AN ( G ( I) 1 /TO ( 1 ) »UT 
T ( 1 ) =UL06 ( TG ( I ) ) /PT 
993 CONTINUE 

.. H I T E ( 6 , 9 9 b ) 

99b FUnMATI* nuM RijOTS S«| 

00996 1=1 , N 

*i< I TE l 6 , 997 I Y ( I I , X ( l I 

997 FORMAT (3rt S» D29.12.iiH + *Nu - J U29.I21 
996 CONTINUE 

CALL «TP0LY(N,B,b0,l.E-2b,HTK,RTI,C0NV,A,o,C,D,E) 
flH 1 TE ( 6 , 99B ) 

998 FUKrlATI' OLNOM ROOTS Z»l 
P 0 9 9 9 1 = 1 , N 

ARITE(6,50UI RTR ( I ) ,RT I l I ) v 

SOO FORMAT I 3H 2« -029.l2.tlM ♦ aNO - J 029.121 

999 CUNTINUE 

U U b u 2 I * 1 > N 

T66II)«OSORT(RTR(Il*KTR(II.RTl(l)«KTlun 
66 1 1 )=RT I II >/R TH 1 1) 

XX I I ) = DA T An (66 ( I I) / TOG ( l I »0 T 
Y Y ( I 1 -DlOG ( TGG ( I > > yOT 
SO 2 CONTINUE 

wHi T E ( 6 , SU3 ) 

SiJ3 TOKnATI' OEMIH ROOTS s'l 
0 U b 0 9 1 ■ 1 , N 

«»R I r e ( 6 , b 0 6 ) Y Y ( I ) , XX I I I 
bU6 F0 RMaT(3h S« 029.12,Iih ♦ A NO - J 029.12) 
bU 9 CONTINUE 
60 TO 1 

bO CaLL PR I NT ( 1 2HENU OF RUN ,1.A) 

STOP 
E N 0 
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